Time Series Modeling

The Challenge

Finance has the responsibility to forecast, plan and allocate resources in order to successfully deliver the best treatments to our patients.

• The first challenge consist of deliver sales forecast for 2018 per cluster, per brand and per month.

• The second challenge is to provide the optimal allocation of their resources for the year 2018.

Sales Forecasting

The first problem deals with generating the forecasts of 75 time series. Each time series corresponds to the demand for a specific product (brand) in a geographic area (cluster). The time series are in monthly granularity.

We have been provided with the demand history for each of the series, so that by training statistical models in the past, we should be able to model behaviors that are likely to be repeated in the future. Therefore, with time series models we should be able to tackle the problem.

There are different techniques to tackle this problem, the “bottom-up” and the “top-down” approach.

A simple method for generating coherent forecasts is the “bottom-up” approach. This approach involves first generating forecasts for each series at the bottom level, and then summing these to produce forecasts for all the higher level series. Other approaches are the “top-down” whose involve first generating forecasts for the total series, and then disaggregating these down to all the levels.

The method chosen for this problem is the “bottom-up” approach because with this method we do not loss any information due to aggregations and because with only 75 series is a feasible approach.

Recall how the data looks like

id	cluster	brand	date	sales_1	sales_2	investment_1	investment_2	investment_3	investment_4	investment_5	investment_6
Cluster 9__Others	Cluster 9	Others	2018 Oct	0	0	22.2	0.0	0.0	-1.8	0	0.0
Cluster 9__Others	Cluster 9	Others	2018 Nov	0	0	-5.5	-7.0	-1.6	-0.6	0	-2.4
Cluster 9__Others	Cluster 9	Others	2018 Dec	0	0	-7.2	-3.1	-1.6	-0.6	0	-2.4

The training process will be divided into the following steps:

1. Define a list of models to use.
2. Train the models using the training data
3. Create the forecasts for the validation data
4. Evaluate the performance of the models

Evaluation metric

The competition metric is the Mean Absolute Percentage Error (MAPE).

\[ MAPE = \frac{100}{h} \sum^{h}_{t=1}\frac{y_{t} - \hat{y}_{t} }{y_{t}} \]

where \(y_t\) is the actual value and \(\hat{y}_t\) is the predicted value.

Percentage errors have the advantage of being unit-free, and so are frequently used to compare forecast performances between data sets.

Measures based on percentage errors have the disadvantage of being infinite or undefined if \(y_t=0\) for any \(t\) in the period of interest, and having extreme values if any \(y_t\) is close to zero. Makridakis (1993) said that “equal errors above the actual value result in a greater APE than those below the actual value”. He provided an example where \(y_t = 150\) and \(\hat{y}_t = 100\), so that the relative error is \(\frac{50}{150} = 0.33\), in contrast to the situation where \(y_t = 100\) and \(\hat{y}_t = 150\), when relative error would be \(\frac{50}{100} = 0.5\). Thus, MAPE puts a heavier penalty on negative errors (when \(y_t\) < \(\hat{y}_t\)) than on positive errors.

In the x-axis we have the absolute error and in the y-axis de MAPE. Absolute error are the number of units above or below the real value, i.e if the real value is 100 and the error is 50 then the predicted value has been 150. For the same error, as lower the true vale it is, the bigger the error we have.

Split data

As we are dealing with monthly time series data we will split the series using time rather than randomly as we usually do when training machine learning models. We will keep apart the last year of the data.

The datasets will be as follows:

[1] "Train data period: 2012 Jan - 2016 Dec"
[1] "Valid data period: 2017 Jan - 2017 Dec"
[1] "Test data period: 2018 Jan - 2018 Dec"

Baseline models

Before doing nothing, let’s compute some statistical baselines for all the series, that is, statistical models that just repeat the last value (NAIVE) or the last period of samples (SNAIVE).

The models that we are going to use as a benchmarks are:

• NAIVE: this model replicate the last value of the sales data. This model is equivalent to an \(ARIMA(0,1,0)\) model.

• SNAIVE: this model replicates the last period of samples from train data. This model is equivalent to an \(ARIMA(0,0,0)(0,1,0)_m\) model, where m is the seasonal period (m = 12 in monthly data.)

# Baseline models
models <- list(
  naive    = NAIVE(sales_2),
  snaive_1 = SNAIVE(sales_2 ~ lag("year") ),
  snaive_2 = SNAIVE(sales_2 ~ lag("year") + drift() )
)

Fit baseline models

fit_baselines <- train_data %>% 
  model(!!!models)

Generate forecast for all the models. Three different forecast will be generated for each serie.

fc_baselines <- fit_baselines %>% 
  forecast(new_data = valid_data)

The metric that we will use to evaluate the performance of each model is the MAPE but I will also include the RMSSE and CRPS for completeness. For both, the RMSE and the MAPE, the smaller the better.

The best model has been the naive with a MAPE of 26.9, followed by the two variants of the snaive. An error of a MAPE of 26.9 means that the model has an error of 26.9%.

.model	RMSSE	MAPE	CRPS
naive	17.1	26.9	1224.8
snaive_1	17.3	29.9	1174.4
snaive_2	17.3	30.1	1197.3

Now that we have the forecast and we know what has been the best model, let’s visualize the forecast to get insight into the data.

Although the snaive model has performed better overall, in some series it performs worse than the snaive.

We also find differences between the snaive models. In the following plot we can see how the snaive with drift is under-forecasting while the other, without drift, is providing an acceptable forecast.

Sometimes no model is good.

Statistical Benchmarks

Exponential smoothing and ARIMA models are the two most widely used approaches to time series forecasting, and provide complementary approaches to the problem. While exponential smoothing models are based on a description of the trend and seasonality in the data, ARIMA models aim to describe the autocorrelations in the data.

• ARIMA: autoregressive integrated moving average model
• ETS: exponential smoothing models are based on a description of the trend and seasonality in the data

When there are long seasonal periods, a dynamic regression with fourier terms is often better than ARIMA and ETS models.

For example, daily data can have annual seasonality of length 365, weekly data has seasonal period of approximately 52, while half-hourly data can have several seasonal periods, the shortest of which is the daily pattern of period 48.

Despite we do not have long seasonal periods here, this models will be included in our list of models.

• FOURIER: Linear regression model with errors modeled with an arima model to remove the autocorrelation.

Additionally to this models, it has also been included a combination forecast of ARIMA and ETS models which has been called combn.

models_stat <- list(
  # Statistical benchmarks
  arima = ARIMA(sales_2),
  ets   = ETS(sales_2),
  
  # Using Fourier terms and ARIMA errors for forecasting
  `K = 1` = ARIMA(sales_2 ~ fourier(K = 1) + PDQ(0, 0, 0)),
  `K = 2` = ARIMA(sales_2 ~ fourier(K = 2) + PDQ(0, 0, 0)),
  `K = 3` = ARIMA(sales_2 ~ fourier(K = 3) + PDQ(0, 0, 0)),
  `K = 4` = ARIMA(sales_2 ~ fourier(K = 4) + PDQ(0, 0, 0)),
  `K = 5` = ARIMA(sales_2 ~ fourier(K = 5) + PDQ(0, 0, 0)),
  `K = 6` = ARIMA(sales_2 ~ fourier(K = 6) + PDQ(0, 0, 0))
)

The best model is the combn, followed by the ARIMA and the K = 2 models. The error has been reduced by 26.4% with respect to the naive model.

.model	RMSSE	MAPE	CRPS
combn	16.7	19.8	797.7
arima	16.9	20.4	827.0
K = 2	16.8	21.6	875.1
ets	16.8	21.8	841.9
K = 5	16.9	21.8	926.6
K = 4	17.0	21.8	935.6

In the example below, neither model has generated a good forecast, in fact, it seems that the ets model only has captured the trend.

In this other case, both forecasts are similar to each other and the snaive_1

Again, none of the models has been able to capture the big drop off of 2017.

External Information

The time series models used so far allow for the inclusion of information from past observations of a series, but not for the inclusion of other information that may also be relevant. For example, the effects of the investments may explain some of the historical variation and may lead to more accurate forecasts

# Using investments as predictors
models_xreg <- list(
  xreg_1 = ARIMA(sales_2 ~ investment_1),
  xreg_2 = ARIMA(sales_2 ~ investment_1 + investment_2),
  xreg_3 = ARIMA(sales_2 ~ investment_1 + investment_2 + investment_3),
  xreg_4 = ARIMA(sales_2 ~ investment_1 + investment_2 + investment_3 + 
                   investment_4 + investment_5 + investment_6)
)

Overall accuracy is not as good as combn, arima orets models trained on statistical benchmarks section.

.model	RMSSE	MAPE	CRPS
xreg_3	17.0	23.9	1065.3
xreg_1	17.1	24.1	958.7
xreg_2	17.1	24.9	1047.7
xreg_4	17.2	25.7	1119.4

Model K = 4 overfit the data and generate a forecast in the opposite direction. In contrast, the K = 2 is able to capture a more general pattern. This is a good example that more complex models do not always indicate better results.

All the forecasts are quite similar in this case

Apparently, the investments neither explain the drop off of 2017.

For cases like this, we are going to introduce a new set of models that will help us to improve this structural changes.

Piecewise regression

We have seen that, in some cases, the trend is not linear and the models are not able to model the series correctly. Piece-wise models divide the series into several segments in such a way that they are able to detect the trend correctly.

Brand analysis

Performing the analysis for each series is a complicated and time-consuming task. For this reason, many times higher levels of aggregation are analyzed and then these conditions or parameters are applied to lower levels series.

The plot of the Brand Group 31 reveals two different periods. There is a clear uptrend up to the end of 2015. After 2015 there is a decrease in the growing speed and the trend slow down a bit. To account for these changes, we specify the date 2015-01 as knots.

# Select a brand a fit a model
current_brand <- "Brand Group 31"

# Fit simple `lm` and `piecewise` model
fit_trends <- by_brand_tsibble %>% 
  filter(brand == current_brand) %>% 
  model(
    lm = TSLM(sales_2 ~ trend()),
    piecewise = TSLM(sales_2 ~ trend(
      knots = tsibble::yearmonth("2015-01")))
  )

In the plot below the relationship between the real and adjusted values by the model is shown. The slashed line represents the perfection, in such a way that the further the points are from the line, the worse the model is. In the lower left corner, the values of the piecewise model are much better adjusted to the reality while those of the lm are farther away. In general, there is a bigger dispersion in the lm model than in the piecewise.

The models are defined manually for each brand.

get_models_ <- function(brand) {
  # Define step-wise models
  mdl <- null_model()
  if (brand == "Brand Group 17") {
    mdl  <- TSLM(sales_2 ~ trend(knots = tsibble::yearmonth("2015-03")) + season())
  } else if (brand == "Brand Group 24") {
    ...
  }
}

In brands such as Brand Group 17 or Brand Group 24 the effect is remarkable. It is also in others, such as the Brand Group 31 as we have just seen in the previous plots.

All series analysis

Now that we have identified the knots of each brand is time to apply this models into the product level series.

fit_stepwise <- map_dfr(
  .x = setNames(brands, brands), 
  .f = ~ train_data %>% 
    filter(brand == .x) %>% 
    model(!!!get_models(.x)["stepwise"]), 
  .id = "brand"
)

Overall accuracy is not as good as combn, arima orets models trained on statistical benchmarks section

.model	RMSSE	MAPE	CRPS
stepwise	16.9	24	1077.7

Let’s look at some of the worst forecast. The arima and ets models will be included as a reference. In that way, we can understand how the piecewise model is modeling the data. In the following case, it seems that the model has failed to capture the trend correctly and is under-forecasting much more than the other two models.

What we can do now? Let’s plot the trend to understand why this is happening. The model fits very well the trend up to 2016 but, after that period, the trend is still in a downtrend and it should remain almost flat. The vertical dashed lines indicate the knots defined within the model.

Finally, we show the forecast for the series that we have been looking in the previous sections.

fit_stepwise %>% 
  filter(id == "Cluster 2__Others") %>%
  left_join(select(fit_statistical, id, ets, arima), by = "id") %>% 
  forecast(new_data = valid_data) %>% 
  plot_forecast(facets = ".model", ncol = NULL, level = NULL, title = "") +
  facet_wrap(~ id)

Selecting best model

Throughout the notebook we have seen how each model is better suited to one specific case, in some cases a model as simple as the naive give us the best results, however, in other cases we need more complex models to capture the signal of the series. So, why not select the best model for each series instead of always using the same model for everything?

# Combine all accuracies together
acc_all <- 
  bind_rows(acc_baselines, acc_statistical, acc_xreg, acc_stepwise) %>% 
  add_rank(var = "MAPE")

# Select best model for each serie
acc_best_models <- acc_all %>% 
  arrange(id, rank) %>% 
  top_n(-rank, n = 1) %>% 
  dplyr::group_by(id) %>% 
  dplyr::slice(1) %>% 
  ungroup()

The error has been reduced by an 47% with respect to the naive and 28% with respect to the combn model.

acc_best_models %>% 
  summarise_at(vars(RMSSE:CRPS), mean) %>% 
  add_column(.model = "Best model", .before = 1) %>% 
  arrange(MAPE) %>% 
  to_html(digits = 1)

.model	RMSSE	MAPE	CRPS
Best model	16.5	14.1	602.1

How many times each model has been selected? stepwise and arima are the two best models followed by xreg_4. Recall that this model uses all the investments provided.

In the same way that we have been doing, we are going to show the best and worst predictions to get an idea of what the predictions are like.

# Combine all forecasts
fc_all_models <- bind_rows(fc_baselines, fc_statistical, fc_xreg, fc_stepwise)

fc_best_models <- fc_all_models %>% 
  semi_join(acc_best_models, by = c(".model", "id"))

Sometimes, any model is able to predict the data accurately. Might be, because the series are white noise or the signal-to-noise ratio is very low. Other times, there are very abrupt changes that cannot be explained by any predictor (past sales, investments, etc.)

Drill down accuracy

Forecasts are often required for all levels of the series, and it is natural that the forecasts to add up in the same way as the data. This is important because the forecast will be analyzed from different points of view and all of them must be aligned. For example, a brand manager might be interested in the forecast of all of their products but it also would might get a view at country level of their brand. At the same time, the CEO of the company, might be interested in higher levels of aggregation.

The easiest brand to forecast is Brand Group 96, 97 and the hardest is Brand Group 24 which is not a surprise given it’s highest growth.

.model	brand	RMSSE	MAPE	CRPS
Best model	Brand Group 96, 97	0.5	6.6	208.9
Best model	Others	0.4	6.8	1315.0
Best model	Brand Group 31	0.3	6.9	765.8
Best model	Brand Group 41	0.5	7.7	470.8
Best model	Brand Group 51, 73, 90	0.4	8.2	281.9
Best model	Brand Group 17	0.5	9.0	362.4
Best model	Brand Group 30	0.5	9.1	444.1
Best model	Brand Group 36	1.2	19.2	178.2
Best model	Brand Group 24	1.4	20.1	227.3

In the same way, the easiest cluster is Cluster 8 and the hardest is Cluster 7

.model	cluster	RMSSE	MAPE	CRPS
Best model	Cluster 8	0.6	6.2	270.8
Best model	Cluster 2	0.5	6.3	512.4
Best model	Cluster 10	0.5	6.6	479.5
Best model	Cluster 9	0.5	8.1	168.1
Best model	Cluster 1	0.3	8.9	386.3
Best model	Cluster 3	0.5	9.7	1018.2
Best model	Cluster 4	0.6	16.0	644.1
Best model	Cluster 5	0.9	16.4	208.9
Best model	Cluster 7	0.9	27.4	275.4

Make forecasts

The last step is generate the forecast for all the 75 series. Now, we must estimate the models again but using all the data available.

# Combine all three list of models
models_list <- c(models, models_stat, models_xreg)

# Fit models using all train data
fit_all_models_1 <- train_data %>%
  bind_rows(valid_data) %>%
  model(!!!models_list) %>%
  mutate(combn = (arima + ets) / 2)

# Train piece-wise models
fit_all_models_2 <- map_dfr(
  .x = setNames(brands, brands),
  .f = ~ train_data %>%
    bind_rows(valid_data) %>%
    filter(brand == .x) %>%
    model(!!!get_models(.x)["stepwise"]),
  .id = "brand"
)

# Combine all models together
fit_all_models <- fit_all_models_1 %>% 
  left_join(fit_all_models_2, by = "id")

# Save models
write_rds(fit_all_models, 
  file.path(path_models, "fit_all_models.rds"), 
  compress = "gz")

Resource Allocation

Resource allocation methods can be used in almost every aspect of supply chain management. The common goal is allocate the resources effectively (i.e., maximizing the revenue, minimizing the cost, or optimizing the utilization sequence) while satisfying certain constraints (e.g., resource availability, customer service level, back order level, delivery window). In this case, resources refers to investments.

In this second part of the challenge, Novartis ask to the participants to provide the optimal allocation of their resources for the year 2018, as well as, showing differences of allocations by cluster/brand/investment type and impact on sales. Investments are now compressed into Investment 1, Investment 2 and Others.

The template file provided by the organization looks as follow:

cluster	brand	function	OptimiseInvestment	ForecastedSales	OptimiseSales
Cluster 1	Brand Group 17	Investment 1	NA	NA	NA
Cluster 1	Brand Group 17	Investment 2	NA	NA	NA
Cluster 1	Brand Group 17	others	NA	NA	NA

Create the new variable investment_others, and remove unnecessary columns.

alloc_df <- all_levels_tbl %>% 
  mutate(
    investment_others = investment_3 + investment_4 + 
      investment_5 + investment_6
  ) %>% 
  select(-c(investment_3, investment_4, investment_5, investment_6))

id	cluster	brand	date	sales_1	sales_2	investment_1	investment_2	investment_others
Cluster 1__Brand Group 17	Cluster 1	Brand Group 17	2012 Jan	0	0	0	0	0
Cluster 1__Brand Group 17	Cluster 1	Brand Group 17	2012 Feb	0	0	0	0	0
Cluster 1__Brand Group 17	Cluster 1	Brand Group 17	2012 Mar	0	0	0	0	0

Forecasting model

We need a forecasting model with investments as predictors. Ideally, the model obtained in the first challenge is the one that we should use to perform the optimization. However, some of the models are univariate and in any case the variable investment_others has been used. Therefore, a linear regression model will be used, which will use the trend, the seasonality and the investments as predictors.

Note: we should look for the best model possible that contains investments as predictors.

models_alloc <- list(
  lm = TSLM(sales_2 ~ trend() + season() + investment_1 + investment_2 +
              investment_others)
)

Then, train the models and generate the forecasts for the next 12 months

# Train models
fit_alloc <- train_data_alloc %>% 
  model(!!!models_alloc)

# Generate forecasts for the test period
fc_alloc <- fit_alloc %>% 
  forecast(new_data = test_data_alloc)

Simulations

When a closed form is not available to estimate the value of a parameter, either because it does not have a closed form formula or because there are too many factor that influence in the final result, the use of simulation techniques are commonly used to observe the behavior of the variable to be analyzed.

The model will be fed with variations of the investments to obtain different forecasts. Sometimes reducing the amount invested on one side to place it in another can lead to better results, and vice versa.

The steps to follow are:

1. Create a matrix with the weights combinations
2. Calculate the new simulated investments
3. Generate the forecast with the new simulated values

Simulation weights

We create a mesh of values that allow us to evaluate different scenarios. The following restrictions have been applied:

• The weights must add up to 1

• In order to diversify, the maximum allocation will be 80% of the capital.

comb_id	comb_desc	investment_1	investment_2	investment_others
1	[0.00, 0.20, 0.80]	0	0.20	0.80
2	[0.00, 0.25, 0.75]	0	0.25	0.75
3	[0.00, 0.30, 0.70]	0	0.30	0.70
4	[0.00, 0.35, 0.65]	0	0.35	0.65
5	[0.00, 0.40, 0.60]	0	0.40	0.60
6	[0.00, 0.45, 0.55]	0	0.45	0.55

A total of 189 combinations per series will be evaluated.

Create scenarios

To generate the multiple scenarios we must multiply the weights generated in the previous section with the values of the investments. To do this, we will help ourselves with the following function:

get_future_scenerarios <- function(future_data, weights_df) {
  
  # Store id data
  id_df <- future_data %>% 
    select(id, date)
  
  # Calculate total invested in each product-month
  total_invested <- future_data %>% 
    transmute(total = rowSums(across(contains("investment"))))
  
  # Repeat `total_invested` column 3 times (one for each investment)
  inv <- replicate(n = 3, total_invested$total) 
  
  # Loop through weights
  inv_cols <- str_subset(names(weights_df), "investment")
  scenarios_list <- vector("list", length = nrow(weights_df))
  for (i in seq(1, nrow(weights_df))) {
    
    # Weights for this iteration
    w <- weights_df[i, inv_cols]
    w <- as.numeric(w)
    
    # New investments
    new_investments <- t(w * t(inv))
    colnames(new_investments) <- inv_cols

    future_scenario <- id_df %>% 
      bind_cols(weights_df[i, "comb_id"]) %>% 
      bind_cols(as_tibble(new_investments))
    
    scenarios_list[[i]] <- future_scenario
  }
  
  # Return a list with all the simulations
  scenarios_list
}

In gray-blue color we see the possible scenarios of the forecast. Each scenario corresponds to a different weight combination. In orange we have the median and in another shade of blue the 80th quantile.

Find optimal values

We will use the 80th quantile as the best combination. We could use the combination of parameters that has generated the greatest benefits, but it would be select the best optimistic scenario and we prefer to be more conservative.

First, we calculate the total amount of each simulation

simulations_agg <- fc %>% 
  as_tibble() %>% 
  group_by(id, comb_id) %>% 
  summarise(across(contains(".mean") | contains("investment"), sum)) %>%
  ungroup() %>% 
  rename(OptimalSales = .mean) %>% 
  arrange(desc(OptimalSales)) %>% 
  left_join(select(weights_df, comb_id, comb_desc), by = "comb_id")

Then, select the combination that meets our requirements: pct_rank <= 0.8

# Find best combination of weights for each `id`
optimal_simulation_id <- simulations_agg %>% 
  group_by(id) %>%
  summarise(
    comb_id = comb_id,
    OptimalSales = OptimalSales,
    pct_rank = percent_rank(OptimalSales)
  ) %>% 
  # Select 80 percentile
  filter(pct_rank <= 0.8) %>% 
  slice_max(pct_rank) %>%
  # In case of tie, choose the first one. This happens
  # in `Cluster 5__Brand Group 30`
  slice_head(n = 1) %>% 
  ungroup()

Finally, prepare the output data

id	function	OptimiseInvestment	ForecastedSales	OptimalSales
Cluster 3__Brand Group 31	investment_1	-5851	427441	442260
Cluster 3__Brand Group 31	investment_2	-4388	427441	442260
Cluster 3__Brand Group 31	investment_others	-19016	427441	442260
Cluster 2__Brand Group 41	investment_1	-1311	354095	412700
Cluster 2__Brand Group 41	investment_2	-13107	354095	412700
Cluster 2__Brand Group 41	investment_others	-11796	354095	412700

According to the simulation results, we could have had a profit of ~ 500k in 2018.

Year	OptimiseInvestment	ForecastedSales	OptimalSales	Profit
2018	-457078	7366401	7885967	519566

Evolution of investments

The evolution of the investments over the years is shown below.It can be seen how, in general terms, the investments for 2018 are aligned with the historical trends and do not different sharply from the recent past. However, there are some investments that have been affected more than the rest.

In general, it is suggested to invest more in others for the 2018.

What’s next

I will leave some ideas that be explored to expand the analysis

• Remove outliers and clean the data.
• Time series clustering based on features such as trend strength, seasonality, spikiness, linearity, autocorrelation, etc.
• Test if ensembles do improve the accuracy. There are many types of ensembles to explore: combination ensembles, weighted average ensembles based on cross validation errors or even more advanced using regularized linear models to name a few.
• Use more than one fold to select the best models, i.e use time series cross validation.
• Directly find the best model in the first challenge so it can be used in the optimization problem.
• Try others methods such as prophet, machine learning or others.
• Account for uncertainty measures, because all is not about accuracy, prediction intervals are also important.
• Try other optimization methods